# 0.0 Getting started ----
## 0.1 Load libraries ----
library(spatialreg)
library(spdep)
library(tidyverse)
library(gt)
library(ggh4x)
library(pscl)
library(wnominate)
library(rio)
## 0.2 Define functions ----
table.format <- function(data, stat){
  sessions <- as.character(119:144)[-c(6, 11, 16)]
  neigh <- c('_adj','_rook','_queen')
  named_sessions <- c()
  for (s in sessions){for (n in neigh){named_sessions<-append(named_sessions, paste0(s,n))}}
  named_sessions <- set_names(rep(c('Adjacent','Rook','Queen'),23),named_sessions)
  
  tab <- data %>% 
    mutate(value = paste0(format(round(est,3), nsmall=3), ' (',format(round(pvalue,3), nsmall=3), ')'),
           dv = factor(dv, levels = c('coord1d','coord2d','irt1d', 'irt2d','coalitionscore', 'defectionrate'), 
                       labels = c('WNOM-1', 'WNOM-2', 'IDEAL-1', 'IDEAL-2', 'Coal. Score', 'Def. Rate'),
                       ordered = TRUE)) %>%
    arrange(dv) %>% 
    select(., -c(est, pvalue)) %>% 
    pivot_wider(names_from = c(session,neighbor), values_from = value) %>%
    gt() %>%
    tab_spanner_delim(delim = '_') %>% 
    cols_label(.list = named_sessions) %>% 
    cols_width(everything() ~ px(90)) %>% 
    text_transform(locations = cells_column_labels(columns= 'dv'), fn = function(x){ paste0("") }) 
  
  if (stat=='g'){
    tab <- tab %>% 
      tab_header(title = "GLOBAL MORAN’S I RESULTS BY SESSION AND NEIGHBOR DEFINITION", 
                 subtitle = 'rep: Appendix tables D1-D4') %>%
      tab_style(locations = cells_title(groups = c('title', 'subtitle')),
                style = list(cell_text(size='x-large', align = 'left')))%>%
      tab_source_note(source_note = 'Note: Table reports Moran’s I test statistics with p-values in parentheses.')
  } else if (stat=='s'){
    tab <- tab %>% tab_header(title = "SPATIAL LAG RESULTS BY SESSION AND NEIGHBOR DEFINITION", 
                              subtitle = 'rep: Appendix tables D5-D8') %>%
      tab_style(locations = cells_title(groups = c('title', 'subtitle')),
                style = list(cell_text(size='x-large', align = 'left')))%>%
      tab_source_note(source_note = 'Note: Table reports spatial lag (ρ) parameter estimates with p-values in parentheses.')
  }
  return(tab)
}

table.format.mb <- function(data, stat){
  sessions <-  c(127,130,135,139,140)
  neigh <- c('_adj','_rook','_queen')
  named_sessions <- c()
  for (s in sessions){for (n in neigh){named_sessions<-append(named_sessions, paste0(s,n))}}
  named_sessions <- set_names(rep(c('Adjacent','Rook','Queen'),5),named_sessions)
  
  tab <- data %>% 
    mutate(value = paste0(format(round(est,3), nsmall=3), ' (',format(round(pvalue,3), nsmall=3), ')'),
           dv = factor(dv, levels = c('coord1dmb','coord2dmb','irt1dmb', 'irt2dmb','coalitionscoremb', 'defectionratemb'), 
                       labels = c('WNOM-1', 'WNOM-2', 'IDEAL-1', 'IDEAL-2', 'Coal. Score', 'Def. Rate'),
                       ordered = TRUE)) %>%
    arrange(session) %>% 
    select(., -c(est, pvalue)) %>% 
    pivot_wider(names_from = c(session,neighbor), values_from = value) %>%
    arrange(dv) %>%
    gt() %>%
    tab_spanner_delim(delim = '_') %>% 
    cols_label(.list = named_sessions) %>% 
    cols_width(everything() ~ px(90)) %>% 
    text_transform(locations = cells_column_labels(columns= 'dv'), fn = function(x){ paste0("") }) %>%
    sub_missing(columns = everything(), rows = everything(), missing_text = "-")
  
  if (stat=='g'){
    tab <- tab %>% 
      tab_header(title = "GLOBAL MORAN’S I RESULTS BY SESSION AND NEIGHBOR DEFINITION, MEMBER BILLS", 
                 subtitle = 'rep: Appendix table E2') %>%
      tab_style(locations = cells_title(groups = c('title', 'subtitle')),
                style = list(cell_text(size='x-large', align = 'left'))) %>%
      tab_source_note(source_note = 'Note: Table reports Moran’s I test statistics with p-values in parentheses.')
  } else if (stat=='s'){
    tab <- tab %>% tab_header(title = "SPATIAL LAG RESULTS BY SESSION AND NEIGHBOR DEFINITION, MEMBER BILLS", 
                              subtitle = 'rep: Appendix table E3') %>%
      tab_style(locations = cells_title(groups = c('title', 'subtitle')),
                style = list(cell_text(size='x-large',align = 'left'))) %>%
      tab_source_note(source_note = 'Note: Table reports spatial lag (ρ) parameter estimates with p-values in parentheses.')
  }
  return(tab)
}
# 0.3 Load environment ----
load('data/0.0 rep_env.RData')
# 1.0 Paper Figures ----
## 1.1 Statistical models ----
local({
  sessions <- (119:144)[-c(6, 11, 16)]
  neighbor_types <- c('adj_def','rook_def','queen_def')
  varlist <- c('coord1d', 'coord2d', 'coalitionscore', 'defectionrate', 'irt1d', 'irt2d')
  Lmorans <<- spatlag <<- Gmorans <<- data.frame(matrix(nrow=0, ncol=5))
  colnames(spatlag) <<- colnames(Gmorans) <<- c('session', 'neighbor', 'dv', 'est', 'pvalue')
  colnames(Lmorans) <<- c('session','neigh','dv', 'pvalue')
  spatial_lag_formulas <- c(' ~ sdp + sdu + ip + wl + pm', 
                            ' ~ sdpp + sdu + ip + wl', 
                            ' ~ sdu + pp + sdpp + ip + sdauf',
                            ' ~ pp + lp + ip + sdauf', ' ~ fl + s + sf + vg',
                            ' ~ f + s + sf + vg', ' ~ hr + s + sf + u + vg', 
                            ' ~ f + p + s + sf + vg')
  summary_stats <<- data.frame(matrix(nrow=0,ncol=14))
  colnames(summary_stats) <<- c('session','obs','wnom','wnom_se','wnom2','wnom2_se','coal',
                                'coal_se','def','def_se',
                                'ideal','ideal_se','ideal2','ideal2_se')
  
  for (i in 1:length(sessions)){
    
    if (i<=5){ neighb_def <- '119.123'
    } else if (i<=9){ neighb_def <- '124.128'
    } else if (i==10){ neighb_def <- '129.130'
    } else if (i==11|i==12){ neighb_def <- '131.132'
    } else if (i==14){neighb_def <- '134.135'
    } else { neighb_def <- as.character(sessions[i])}
    
    adj_def <- eval(as.name(paste0('neighbor.adj_',neighb_def)))
    rook_def  <- eval(as.name(paste0('neighbor.rook_',neighb_def)))
    queen_def <- eval(as.name(paste0('neighbor.queen_',neighb_def)))
    if (sessions[i]<121) { sl_formula <- spatial_lag_formulas[1]
    } else if (sessions[i]<123){ sl_formula <- spatial_lag_formulas[2]
    } else if (sessions[i]==123){ sl_formula <- spatial_lag_formulas[3]
    } else if (sessions[i]==125){ sl_formula <- spatial_lag_formulas[4]
    } else if (sessions[i]<130){ sl_formula <- spatial_lag_formulas[5]
    } else if (sessions[i]<138){ sl_formula <- spatial_lag_formulas[6]
    } else if (sessions[i]<142){ sl_formula <- spatial_lag_formulas[7]
    } else { sl_formula <- spatial_lag_formulas[8] }
    
    for (t in 1:length(neighbor_types)){
      
      if ((sessions[i]==137 | sessions[i]==139) & 
          neighbor_types[t] == 'adj_def'){
        session_data <- paste0('session',sessions[i],'_adj')
      } else if ((sessions[i]==137 | sessions[i]==139) &
                 neighbor_types[t] != 'adj_def'){
        session_data <- paste0('session',sessions[i],'_QueenRook')
      } else {session_data <- paste0('session',sessions[i])}
      
      if (t==2) {
        summary_stats[i,1] <<- sessions[i]
        summary_stats[i,2] <<- nrow(eval(as.name(session_data)))
      } else {NULL}
      
      for (v in 1:length(varlist)){
        moran_stat <- ifelse(i==1 & varlist[v]==varlist[3],'cscore', varlist[v])
        ### 1.1 Global Moran's ----
        temp_globalmorans <- data.frame(matrix(nrow=1, ncol=5))
        colnames(temp_globalmorans) <- c('session', 'neighbor', 'dv', 'est', 'pvalue')
        temp_globalmorans$session <- sessions[i]
        temp_globalmorans$neighbor <- sub('_def', '', neighbor_types[t])
        temp_globalmorans$dv <- varlist[v]
        temp_gmoran_results <- moran.test(
          eval(as.name(session_data))[[moran_stat]],
          eval(as.name(neighbor_types[t])), na.action=na.omit,
          randomisation=TRUE, zero.policy=TRUE)
        temp_globalmorans$est <- temp_gmoran_results$est[[1]]
        temp_globalmorans$pvalue <- temp_gmoran_results$p.value
        Gmorans <<- rbind(Gmorans, temp_globalmorans)
        ## 1.1 Local Moran's ----
        temp_localmorans <- data.frame(matrix(nrow=nrow(eval(as.name(session_data))), ncol=4))
        colnames(temp_localmorans) <- c('session','neigh','dv','pvalue')
        temp_localmorans$pvalue <- localmoran(
          eval(as.name(session_data))[[moran_stat]],
          eval(as.name(neighbor_types[t])),
          zero.policy = TRUE, na.action = na.exclude, 
          conditional = FALSE, alternative = 'greater')[,5]
        temp_localmorans$session <- sessions[i]
        temp_localmorans$neigh <- sub('_def', '', neighbor_types[t])
        temp_localmorans$dv <- varlist[v]
        Lmorans <<- rbind(Lmorans, temp_localmorans)
        ## 1.1 Spatial lag ----
        temp_spatlag <- data.frame(matrix(nrow=1, ncol=5))
        colnames(temp_spatlag) <- c('session', 'neighbor', 'dv', 'est', 'pvalue')
        temp_spatlag$session <- sessions[i]
        temp_spatlag$neighbor <- sub('_def', '', neighbor_types[t])
        temp_spatlag$dv <- varlist[v]
        temp_spatlag_results <- summary(lagsarlm(
          eval(as.formula(paste0(moran_stat, sl_formula))),
          data = eval(as.name(session_data)), 
          listw = eval(as.name(neighbor_types[t])),
          na.action = na.omit, zero.policy=TRUE))
        temp_spatlag$est <- temp_spatlag_results$rho[[1]]
        temp_spatlag$pvalue <- temp_spatlag_results$LR1$p.value[[1]]
        spatlag <<- rbind(spatlag, temp_spatlag)
        ## 1.1 Summary statistics ----
        if (t==2) {
          summary_stats[i,(v*2)+1] <<- mean(eval(as.name(session_data))[moran_stat][[1]], na.rm=TRUE)
          summary_stats[i,(v*2)+2] <<- sd(eval(as.name(session_data))[moran_stat][[1]], na.rm=TRUE)
        } else {NULL}
      }
    }
  }
})
## 1.2 Number of bills and motions (Figure 2) ----
import('data/0.1 malaskr.dta') %>% 
  filter(session >=119  & (type=='l' | type=='a')) %>%
  count(type, session) %>% mutate(type=ifelse(type=='l', 'bills','motions')) %>%
  ggplot(.) + geom_col(aes(session,n,fill=type), position='dodge') + 
  scale_x_continuous(breaks = seq(119,144,2)) + 
  scale_y_continuous(limits = c(0,325), breaks = seq(0,300,100)) +
  labs(title = 'NUMBER OF BILLS & MOTIONS PER SESSION', 
       subtitle = 'ref: Figure 2')
## 1.3 Number of votes per session (Figure 3) ----
local({
  sessions <- 119:144
  vps <<- data.frame(matrix(nrow=length(sessions), ncol=2))
  colnames(vps) <<- c('session','num_votes')
  for (i in 1:length(sessions)) {
    dat <- read.csv(paste0("data/0.2 atkvskr_",sessions[i],".csv"))
    n <- length(unique(dat$v1))
    vps[i,1] <<- sessions[i]
    vps[i,2] <<- n
  }
})
ggplot(vps) + geom_col(aes(session,num_votes)) + scale_x_continuous(breaks = seq(119,144,2)) +
  labs(title = 'NUMBER OF VOTES PER SESSION', 
       subtitle = 'ref: Figure 3')
## 1.4 Global Moran's figure (Figure 5) ----
Gmorans %>%
  mutate(dv = factor(dv, levels = c('coord1d','irt1d','coalitionscore', 
                                    'coord2d','irt2d', 'defectionrate'), 
                     labels = c('WNOM-1', 'IDEAL-1', 'Coalition Score', 
                                'WNOM-2', 'IDEAL-2', 'Defection Rate')),
         neighbor = factor(neighbor, levels = c('adj','rook','queen'), 
                           labels = c('Adjacent','Rook','Queen'))) %>%
  ggplot(., aes(pvalue, reorder(session, -session), shape = neighbor)) +
  geom_point() + 
  labs(x = "p-value", y = "Session", shape="Neighbor Definition: ", 
       title = 'GLOBAL MORAN’S I SIGNIFICANCE LEVELS',
       subtitle = 'ref: Paper Figure 5') +
  geom_vline(xintercept=0.05,linetype = "dashed") + 
  facet_wrap(~ dv) + theme(legend.position = "bottom")
# Statistics referenced in narrative (Global Moran's)
sum(!is.na(Gmorans$pvalue))-sum(Gmorans$pvalue<0.05, na.rm = TRUE)
Gmorans %>% group_by(session) %>% summarise(n=sum(pvalue<0.05)) %>% filter(n==0) %>% count()
Gmorans %>% group_by(dv) %>% summarise(n=sum(pvalue<0.05))
Gmorans %>% group_by(neighbor) %>% summarise(n=sum(pvalue<0.05))
#ggsave("figs/GlobalMorans.pdf", plot=last_plot(), width=8.5, height=11)
## 1.5 Local Moran's figure (Figure 6) ----
Lmorans %>% 
  filter(session==119) %>%
  mutate(dv = factor(dv, levels = c('coord1d','coord2d','irt1d', 'irt2d',
                                    'coalitionscore', 'defectionrate'), 
                     labels = c('WNOM-1', 'WNOM-2', 'IDEAL-1', 'IDEAL-2', 
                                'Coal. Score', 'Def. Rate'),
                     ordered = TRUE),
         neigh = factor(neigh, levels = c('adj', 'rook', 'queen'), 
                        labels = c('Adjacent', 'Rook', 'Queen'))) %>%
  ggplot(., aes(x=pvalue)) + 
  geom_histogram(color="grey30", fill = "white", binwidth = 0.05, 
                 boundary = 0, closed = "left") + 
  labs(x = "p-value", y = "Count", 
       title = 'LOCAL MORAN’S I SIGNIFICANCE LEVELS FOR THE 119TH SESSION',
       subtitle = 'ref: Paper Figure 6') + 
  geom_vline(xintercept = 0.05, linetype = "dashed", alpha = .5) + 
  facet_nested(neigh ~ dv) + theme_bw() + 
  scale_x_continuous(breaks = c(0.0, 0.5, 1.0)) + 
  scale_y_continuous(breaks = c(0, 5, 10, 15))
# Statistics referenced in narrative (Local Moran's)
sum(!is.na(Lmorans$pvalue))
sum(Lmorans$pvalue<0.05, na.rm = TRUE)
#ggsave("figs/LocalMorans119.pdf", plot=last_plot(), width=11, height=8.5)
## 1.6 Spatial lag figure (Figure 7) ----
spatlag %>% 
  mutate(dv = factor(dv, levels = c('coord1d','irt1d','coalitionscore', 
                                    'coord2d','irt2d', 'defectionrate'), 
                     labels = c('WNOM-1', 'IDEAL-1', 'Coalition Score', 
                                'WNOM-2', 'IDEAL-2', 'Defection Rate')),
         neighbor = factor(neighbor, levels = c('adj','rook','queen'), 
                           labels = c('Adjacent','Rook','Queen'))) %>%
  ggplot(., aes(pvalue, reorder(session, -session), 
                                         shape = neighbor)) +
  geom_point() + 
  labs(x = "p-value", y = "Session", shape="Neighbor Definition: ", 
       title = 'SPATIAL LAG SIGNIFICANCE LEVELS', 
       subtitle = 'ref: Paper Figure 7') +
  geom_vline(xintercept = 0.05, linetype = "dashed") + 
  facet_wrap(~ dv) + theme(legend.position = "bottom")
# Statistics referenced in narrative (Spatial lag)
sum(!is.na(spatlag$pvalue))-sum(spatlag$pvalue<0.05, na.rm = TRUE)
#ggsave("figs/SpatialLag.pdf", plot=last_plot(), width=8.5, height=11)
# 2.0 Appendix Tables/Figures ----
## 2.1 Member Bill statistical models ----
load('data/0.3 member_bills.RData')
local({
  sessions <- c(127,130,135,139,140)
  neighbor_types <- c('adj_def','rook_def','queen_def')
  varlist <- c('coord1dmb', 'coord2dmb', 'coalitionscoremb', 'defectionratemb', 'irt1dmb', 'irt2dmb')
  Lmorans.mb <<- spatlag.mb <<- Gmorans.mb <<- data.frame(matrix(nrow=0, ncol=5))
  colnames(spatlag.mb) <<- colnames(Gmorans.mb) <<- c('session', 'neighbor', 'dv', 'est', 'pvalue')
  colnames(Lmorans.mb) <<- c('session','neigh','dv', 'pvalue')
  spatial_lag_formulas <- c(' ~ fl + s + sf + vg',' ~ f + s + sf + vg', ' ~ hr + s + sf + u + vg')
  
  for (i in 1:length(sessions)){
    
    if (i==1){ neighb_def <- '124.128'
    } else if (i==2){ neighb_def <- '129.130'
    } else if (i==3){ neighb_def <- '134.135'
    } else { neighb_def <- as.character(sessions[i])}
    
    adj_def <- eval(as.name(paste0('neighbor.adj_',neighb_def)))
    rook_def  <- eval(as.name(paste0('neighbor.rook_',neighb_def)))
    queen_def <- eval(as.name(paste0('neighbor.queen_',neighb_def)))
    if (sessions[i]<130){ sl_formula <- spatial_lag_formulas[1]
    } else if (sessions[i]<138){ sl_formula <- spatial_lag_formulas[2]
    } else { sl_formula <- spatial_lag_formulas[3] }
    
    for (t in 1:length(neighbor_types)){
      
      if ((sessions[i]==139) & neighbor_types[t] == 'adj_def'){
        session_data <- paste0('session',sessions[i],'_adjmb')
      } else if ((sessions[i]==139) & neighbor_types[t] != 'adj_def'){
        session_data <- paste0('session',sessions[i],'_QueenRookmb')
      } else {session_data <- paste0('session',sessions[i],'mb')}
      for (v in 1:length(varlist)){
        moran_stat <- varlist[v]
        
        if (sessions[i]==130){ NULL
        } else if (sessions[i]==135 & v %in% c(1:3)){ next
        } else if (v %in% 1:2){ next 
        } else (NULL) 

        ### 2.1 Global Moran's ----
        temp_globalmorans <- data.frame(matrix(nrow=1, ncol=5))
        colnames(temp_globalmorans) <- c('session', 'neighbor', 'dv', 'est', 'pvalue')
        temp_globalmorans$session <- sessions[i]
        temp_globalmorans$neighbor <- sub('_def', '', neighbor_types[t])
        temp_globalmorans$dv <- varlist[v]
        temp_gmoran_results <- moran.test(
          eval(as.name(session_data))[[moran_stat]],
          eval(as.name(neighbor_types[t])), na.action=na.omit,
          randomisation=TRUE, zero.policy=TRUE)
        temp_globalmorans$est <- temp_gmoran_results$est[[1]]
        temp_globalmorans$pvalue <- temp_gmoran_results$p.value
        Gmorans.mb <<- rbind(Gmorans.mb, temp_globalmorans)
        ## 2.1 Local Moran's ----
        temp_localmorans <- data.frame(matrix(nrow=nrow(eval(as.name(session_data))), ncol=4))
        colnames(temp_localmorans) <- c('session','neigh','dv','pvalue')
        temp_localmorans$pvalue <- localmoran(
          eval(as.name(session_data))[[moran_stat]],
          eval(as.name(neighbor_types[t])),
          zero.policy = TRUE, na.action = na.exclude, 
          conditional = FALSE, alternative = 'greater')[,5]
        temp_localmorans$session <- sessions[i]
        temp_localmorans$neigh <- sub('_def', '', neighbor_types[t])
        temp_localmorans$dv <- varlist[v]
        Lmorans.mb <<- rbind(Lmorans.mb, temp_localmorans)
        ## 2.1 Spatial lag ----
        temp_spatlag <- data.frame(matrix(nrow=1, ncol=5))
        colnames(temp_spatlag) <- c('session', 'neighbor', 'dv', 'est', 'pvalue')
        temp_spatlag$session <- sessions[i]
        temp_spatlag$neighbor <- sub('_def', '', neighbor_types[t])
        temp_spatlag$dv <- varlist[v]
        temp_spatlag_results <- summary(lagsarlm(
          eval(as.formula(paste0(moran_stat, sl_formula))),
          data = eval(as.name(session_data)), 
          listw = eval(as.name(neighbor_types[t])),
          na.action = na.omit, zero.policy=TRUE))
        temp_spatlag$est <- temp_spatlag_results$rho[[1]]
        temp_spatlag$pvalue <- temp_spatlag_results$LR1$p.value[[1]]
        spatlag.mb <<- rbind(spatlag.mb, temp_spatlag)
      }
    }
  }
})
## 2.2 Representative scaling: 136th session of AlÞingi (Figure A1) ----
### 2.2a W-NOMINATE ----
local({
  load('data/0.4 scaling_136.RData')
  rc <<- rollcall(data=althingi136[,4:ncol(althingi136)], yea=1, nay=6, missing=c(8,9),
                 notInLegis=NA, legis.names=althingi136$NumericID,
                 vote.names=colnames(althingi136[4:ncol(althingi136)]),
                 legis.data=althingi136[,1:3][-2], vote.data=NULL,
                 desc="136th Althingi Roll Call Data")
  result <- wnominate(rc, ubeta=15, uweights=0.5, dims=2, minvotes=20,
                      lop=0.025, trials=3, polarity=c(1,7), verbose=FALSE)
  WEIGHT <- (result$weights[2])/(result$weights[1])
  X1 <- result$legislators$coord1D
  X2 <- (result$legislators$coord2D)*WEIGHT
  party <<- result$legislators$PartyID
  dev.off()
  #pdf(file = "figs/A1nominate.pdf", width = 7, height = 7)
  par(mar=c(5.1,5.1,4.1,2.1),mgp = c(4, 1, 0))
  plot(X1, X2, main="The 136th Althingi: W-NOMINATE",
       xlab="First Dimension (Liberal - Conservative)
    S = Conservatives, F = Progressives,
    U = SDA, L = Liberal, G = Greens",
       ylab="Second Dimension", xlim=c(-1,1), ylim=c(-1,1), asp=1, type="n")
  points(X1[party == 111], X2[party == 111], pch="S", col="gray33", font=2)
  points(X1[party == 105], X2[party == 105], pch="F", col="gray50", font=2)
  points(X1[party == 112], X2[party == 112], pch="U", col="gray33", font=2)
  points(X1[party == 115], X2[party == 115], pch="G", col="gray50", font=2)
  points(X1[party == 106], X2[party == 106], pch="L", col="gray50", font=2)
  p <- recordPlot()
  dev.off()
  return(p)
})
### 2.2b Bayesian IRT/IDEAL ----
local({
  set.seed(123)
  ideal.unidentified <- ideal(rc, d=2, maxiter=55000, thin=50, burnin=5000, normalize=FALSE, store.item=TRUE)
  ideal.identified <- postProcess(ideal.unidentified, 
                                  constraints=list("395"=c(-0.977, -0.0947), "137"=c(0.1432, 0.9897),"294"=c(0.8797, -0.4756)))
  ideal1<-ideal.identified$xbar[,1]
  ideal2<-ideal.identified$xbar[,2]
  dev.off()
  #pdf(file = "figs/A1irt.pdf", width = 7, height = 7)
  par(mar=c(5.1,5.1,4.1,2.1),mgp = c(4, 1, 0))
  plot(ideal1, ideal2, main="The 136th Althingi: IDEAL",
       xlab="First Dimension (Liberal - Conservative)
    S = Conservatives, F = Progressives,
    U = SDA, L = Liberal, G = Greens",
       ylab="Second Dimension", xlim=c(-2.0,2.0), ylim=c(-2.5,2.5), asp=1, type="n")
  points(ideal1[party == 111], ideal2[party == 111], pch="S", col="gray33", font=2)
  points(ideal1[party == 105], ideal2[party == 105], pch="F", col="gray50", font=2)
  points(ideal1[party == 112], ideal2[party == 112], pch="U", col="gray33", font=2)
  points(ideal1[party == 115], ideal2[party == 115], pch="G", col="gray50", font=2)
  points(ideal1[party == 106], ideal2[party == 106], pch="L", col="gray50", font=2)
  p <- recordPlot()
  dev.off()
  return(p)
})
## 2.3 Summary statistics for vote scores & ideal point estimates (Table A1) ----
local({
  new_row <- data.frame(matrix(nrow=1,ncol=14))
  colnames(new_row) <- colnames(summary_stats)
  new_row[1,1] <-'Total'
  new_row[1,2] <-sum(summary_stats[2])
  for (i in 3:14){new_row[i] <- sum(summary_stats[i][[1]]*(summary_stats[2][[1]])/sum(summary_stats[2]))}
  summary_stats <- rbind(summary_stats,new_row)
  summary_stats[c(3:6,11:14)] <- format(round(summary_stats[c(3:6,11:14)],3),nsmall=3)
  summary_stats[7:10] <- format(round(summary_stats[7:10],2),nsmall=2)
  
  p.sumstats <- summary_stats %>% 
    mutate(WNOM1 = paste0(wnom," (",wnom_se,")"),
           WNOM2 = paste0(wnom2," (",wnom2_se,")"),
           IDEAL1 = paste0(ideal," (",ideal_se,")"),
           IDEAL2 = paste0(ideal2," (",ideal2_se,")"),
           CoalScore = paste0(coal," (",coal_se,")"),
           DefRate = paste0(def," (",def_se,")"))
  p.sumstats <- p.sumstats[-c(3:14)]
  tab <- p.sumstats %>% gt() %>% opt_row_striping() %>% opt_vertical_padding(scale=.1) %>%
    tab_header(title = "SUMMARY STATISTICS FOR VOTE SCORES & IDEAL POINT ESTIMATES", 
              subtitle = 'rep: Appendix table A1') %>%
    tab_style(locations = cells_title(groups = c('title', 'subtitle')),
              style = list(cell_text(size='x-large', align = 'left'))) %>%
    tab_source_note(source_note = 'Note: Means (standard deviations) are reported; observations refer to the number of legislators who cast votes.')
  return(tab)
})
## 2.4 Global Moran's table (Tables D1-D4) ----
table.format(Gmorans, 'g')
## 2.5 Spatial lag table (Tables D5-D8) ----
table.format(spatlag, 's')
## 2.6 Local morans figures (Figure D1) ---- 
local({
  sessions <- (119:144)[-c(6, 11, 16)]
  for (i in 1:ceiling(length(sessions)/4)){
    lb <- i+3*(i-1)
    ub <- min(i+(3*i),length(sessions))
    temp.p <- Lmorans %>% 
      filter(session %in% sessions[lb:ub]) %>%
      mutate(dv = factor(dv, levels = c('coord1d','coord2d','irt1d', 'irt2d','coalitionscore', 'defectionrate'), 
                         labels = c('WNOM-1', 'WNOM-2', 'IDEAL-1', 'IDEAL-2', 'Coal. Score', 'Def. Rate'),
                         ordered = TRUE), 
             neigh = factor(neigh, levels = c('adj', 'rook', 'queen'), 
                            labels = c('Adj.','Rook','Queen'))) %>%
      ggplot(., aes(x=pvalue)) + 
      geom_histogram(color="grey30", fill = "white", binwidth = 0.05, 
                     boundary = 0, closed = "left") + 
      labs(x = "p-value", y = "Count") + 
      geom_vline(xintercept = 0.05, linetype = "dashed", alpha = .5) + 
      facet_nested(session + neigh ~ dv) + theme_bw() + 
      scale_x_continuous(breaks = seq(0,1,.5)) + 
      scale_y_continuous(limits = c(0,30),breaks = seq(0,30,10)) + 
      labs(title = paste0('LOCAL MORAN’S I SIGNIFICANCE LEVELS, ' ,sessions[lb], '-',sessions[ub]))
    assign(paste0('figD1_',i), temp.p, inherits = TRUE)
  }
})
### 2.6a Figure D1 on page S18 ----
figD1_1
#ggsave("figs/LocalMoransD1.pdf", plot=last_plot(), width=8.5, height=11)
### 2.6b Figure D1 on page S19 ----
figD1_2
#ggsave("figs/LocalMoransD2.pdf", plot=last_plot(), width=8.5, height=11)
### 2.6c Figure D1 on page S20 ----
figD1_3
#ggsave("figs/LocalMoransD3.pdf", plot=last_plot(), width=8.5, height=11)
### 2.6d Figure D1 on page S21 ----
figD1_4
#ggsave("figs/LocalMoransD4.pdf", plot=last_plot(), width=8.5, height=11)
### 2.6e Figure D1 on page S22 ----
figD1_5
#ggsave("figs/LocalMoransD5.pdf", plot=last_plot(), width=8.5, height=11)
### 2.6f Figure D1 on page S23 ----
figD1_6
#ggsave("figs/LocalMoransD6.pdf", plot=last_plot(), width=8.5, height=8.25)
## 2.7 Votes on member bills (Table E1) ----
merged <- import('data/0.5 votes_merged.dta')
local({
  # Column 1
  dat <- import('data/0.1 malaskr.dta')
  df <- data.frame(table(dat$session[dat$session>=119 & dat$proposaltype=="Þ"]))
  colnames(df) <- c('session','Member Bills')
  # Column 2
  temp <- merged %>% 
    mutate(session = factor(session, levels = as.character(119:144))) %>% 
    filter(proposaltype=="Þ", pctyea<.9) %>% count(session,.drop = FALSE, name = 'Votes <90% Support')
  df <- merge(df,temp)
  # Column 3
  sessions = c(127,130,135,139,140)
  temp <- merged %>% 
    filter(proposaltype=="Þ", pctyea<1, session %in% as.character(sessions)) %>% 
    count(session,.drop = FALSE, name = 'Votes <100% Support')
  df <- merge(df,temp, all=TRUE)
  # Column 4
  temp <- data.frame(matrix(nrow=length(sessions),ncol=2))
  colnames(temp) <- c('session','Party Votes')
  for (s in 1:length(sessions)){
    dat <- import(paste0('data/0.6 wide_member_bills_',sessions[s],'.dta'))
    vote_cols <- grep("vote_", colnames(dat))
    dat2 <- dat %>% 
      mutate_if(is.numeric, ~na_if(., 8)) %>%
      mutate_if(is.numeric, ~na_if(., 9)) %>% 
      mutate_at(.vars = vote_cols, ~(PartyID*10+.))
    
    if (sessions[s]<=135){ party_codes <- c(101,102,105,106,108,109,111:116)
    } else { party_codes <- c(101,102,104:109,111:116) }
    for (i in 1:length(party_codes)){
      Yvec <- c()
      Nvec <- c()
      for (c in 1:length(vote_cols)){
        Yvec <- append(Yvec, sum(dat2[vote_cols[c]]==as.numeric(paste0(party_codes[i],1)), na.rm = TRUE))
        Nvec <- append(Nvec, sum(dat2[vote_cols[c]]==as.numeric(paste0(party_codes[i],6)), na.rm = TRUE))
      }
      assign(paste0('yea',party_codes[i]),Yvec)
      assign(paste0('nay',party_codes[i]),Nvec)
      assign(paste0('pct',party_codes[i]),Yvec/(Yvec+Nvec))
    }
    if (sessions[s]<=135){ 
      pct_coalition <- (yea105+yea111)/(yea105+nay105+yea111+nay111)
      pct_opposition <- (yea101+yea102+yea106+yea108+yea109+yea112+yea113+yea114+yea115+yea116)/(yea101+yea102+yea106+yea108+yea109+yea112+yea113+yea114+yea115+yea116+nay101+nay102+nay106+nay108+nay109+nay112+nay113+nay114+nay115+nay116)
    } else {
      pct_coalition <- (yea112+yea115)/(yea112+nay112+yea115+nay115)
      pct_opposition <- (yea105+yea107+yea111+yea114)/(yea105+yea107+yea111+yea114+nay105+nay107+nay111+nay114)
    }
    temp[s,1] <- sessions[s]
    temp[s,2] <- sum((pct_coalition>0.5 & pct_coalition<=1.0 & pct_opposition<0.5 & pct_opposition>=0.0)|(pct_opposition>0.5 & pct_opposition<=1.0 & pct_coalition<0.5 & pct_coalition>=0.0))
  }
  df <- merge(df,temp, all=TRUE)
  # Column 5  
  temp <- data.frame(matrix(nrow=length(sessions),ncol=2))
  colnames(temp) <- c('session','Memb. w/ Non-zero Def.Rate')
  for (s in 1:length(sessions)){
    dat <- import(paste0('data/0.6 wide_member_bills_',sessions[s],'.dta'))
    vote_cols <- grep("vote_", colnames(dat))
    dat <- dat %>% 
      mutate_if(is.numeric, ~na_if(., 8)) %>%
      mutate_if(is.numeric, ~na_if(., 9))
    
    if (sessions[s]<=135){ party_codes <- c(101,102,105,106,108,109,111:116)
    } else { party_codes <- c(101,102,104:109,111:116) }
    
    defections <- c()
    for (i in 1:length(party_codes)){
      dat2 <- dat %>% filter(PartyID==party_codes[i])
      
      party_vec <- c()
      for (c in 1:length(vote_cols)){
        pct <- sum(grepl('1$',dat2[vote_cols[c]][,1]))/sum(grepl('1$',dat2[vote_cols[c]][,1]),grepl('6$',dat2[vote_cols[c]][,1]))
        maj <- case_when(pct>.5 ~ 1,
                         pct<.5 ~ 6)
        party_vec <- append(party_vec, maj)
      }
      for (r in 1:nrow(dat2)){
        vote_vector <- as.numeric(dat2[r,4:length(dat2)]) 
        n_defections <- sum( as.numeric(dat2[r,4:length(dat2)]) != party_vec, na.rm = TRUE )
        if (length(vote_vector)==sum(is.na(vote_vector))){n_defections <- -1}
        defections <- append(defections, n_defections)
      }
    }
    temp[s,1] <- sessions[s]
    temp[s,2] <- sum(defections>0)
  }
  df <- merge(df,temp, all=TRUE)
  # Table
  tab <- df %>% gt()%>% opt_vertical_padding(scale=.1) %>% 
    sub_missing(columns = everything(), rows = everything(), missing_text = "") %>% 
    tab_header(title = "VOTES ON MEMBER BILLS", subtitle = 'rep: Appendix table E1') %>%
    tab_style(locations = cells_title(groups = c('title', 'subtitle')),
              style = list(cell_text(size='x-large', align = 'left'))) %>%
    cols_width(everything() ~ px(75)) %>% opt_row_striping()
  return(tab)
})
## 2.8 Global Moran's I statistics, Member bills (Figure E1) ----
Gmorans.mb %>%
  mutate(dv = factor(dv, levels = c('coord1dmb','irt1dmb','coalitionscoremb', 
                                    'coord2dmb','irt2dmb', 'defectionratemb'), 
                     labels = c('WNOM-1', 'IDEAL-1', 'Coalition Score', 
                                'WNOM-2', 'IDEAL-2', 'Defection Rate')),
         neighbor = factor(neighbor, levels = c('adj','rook','queen'), 
                           labels = c('Adjacent','Rook','Queen'))) %>%
  ggplot(., aes(pvalue, reorder(session, -session), shape = neighbor)) +
  geom_point() + 
  labs(x = "p-value", y = "Session", shape="Neighbor Definition: ", 
       title = 'GLOBAL MORAN’S I SIGNIFICANCE LEVELS, MEMBER BILLS',
       subtitle = 'ref: Appendix Figure E1') +
  geom_vline(xintercept=0.05,linetype = "dashed") + 
  facet_wrap(~ dv) + theme(legend.position = "bottom")
#ggsave("figs/GlobalMoransMBE1.pdf", plot=last_plot(), width=948, height=601, units = 'px', scale = 3)
## 2.9 Spatial lag results, Member bills (Figure E2) ----
spatlag.mb %>% 
  mutate(dv = factor(dv, levels = c('coord1dmb','irt1dmb','coalitionscoremb', 
                                    'coord2dmb','irt2dmb', 'defectionratemb'),  
                     labels = c('WNOM-1', 'IDEAL-1', 'Coalition Score', 
                                'WNOM-2', 'IDEAL-2', 'Defection Rate')),
         neighbor = factor(neighbor, levels = c('adj','rook','queen'), 
                           labels = c('Adjacent','Rook','Queen'))) %>%
  ggplot(., aes(pvalue, reorder(session, -session), 
                shape = neighbor)) +
  geom_point() + 
  labs(x = "p-value", y = "Session", shape="Neighbor Definition: ", 
       title = 'SPATIAL LAG SIGNIFICANCE LEVELS, MEMBER BILLS', 
       subtitle = 'ref: Appendix Figure E2') +
  geom_vline(xintercept = 0.05, linetype = "dashed") + 
  facet_wrap(~ dv) + theme(legend.position = "bottom")
#ggsave("figs/SpatialLagMBE2.pdf", plot=last_plot(), width=948, height=601, units = 'px', scale = 3)
## 2.10 Local Moran's I statistics, Member bills (Figure E3) ----
Lmorans.mb %>% 
  filter(session %in% 127:139) %>%
  mutate(dv = factor(dv, levels = c('coord1dmb','coord2dmb','irt1dmb', 'irt2dmb','coalitionscoremb', 'defectionratemb'), 
                     labels = c('WNOM-1', 'WNOM-2', 'IDEAL-1', 'IDEAL-2', 'Coal. Score', 'Def. Rate'),
                     ordered = TRUE), 
         neigh = factor(neigh, levels = c('adj', 'rook', 'queen'), 
                        labels = c('Adj.','Rook','Queen'))) %>%
  ggplot(., aes(x=pvalue)) + 
  geom_histogram(color="grey30", fill = "white", binwidth = 0.05, 
                 boundary = 0, closed = "left") + 
  labs(x = "p-value", y = "Count") + 
  geom_vline(xintercept = 0.05, linetype = "dashed", alpha = .5) + 
  facet_nested(session + neigh ~ dv) + theme_bw() + 
  scale_x_continuous(breaks = c(0.0, 0.5, 1.0)) + 
  scale_y_continuous(limits = c(0,40), breaks = seq(0,40,10)) + 
  labs(title = 'LOCAL MORAN’S I SIGNIFICANCE LEVELS, 127-140', 
       subtitle = 'Ref: Appendix Figure E3, page S28')
#ggsave("figs/LocalMoransMBE3_1.pdf", plot=last_plot(), width=8.5, height=11)
Lmorans.mb %>% 
  filter(session==140) %>%
  mutate(dv = factor(dv, levels = c('coord1dmb','coord2dmb','irt1dmb', 'irt2dmb','coalitionscoremb', 'defectionratemb'), 
                     labels = c('WNOM-1', 'WNOM-2', 'IDEAL-1', 'IDEAL-2', 'Coal. Score', 'Def. Rate'),
                     ordered = TRUE), 
         neigh = factor(neigh, levels = c('adj', 'rook', 'queen'), 
                        labels = c('Adj.','Rook','Queen'))) %>%
  ggplot(., aes(x=pvalue)) + 
  geom_histogram(color="grey30", fill = "white", binwidth = 0.05, 
                 boundary = 0, closed = "left") + 
  labs(x = "p-value", y = "Count") + 
  geom_vline(xintercept = 0.05, linetype = "dashed", alpha = .5) + 
  facet_nested(session + neigh ~ dv) + theme_bw() + 
  scale_x_continuous(breaks = c(0.0, 0.5, 1.0)) + 
  scale_y_continuous(limits = c(0,40), breaks = seq(0,40,10)) + 
  labs(title = paste0('LOCAL MORAN’S I SIGNIFICANCE LEVELS, 127-140'),
       subtitle = 'Ref: Appendix Figure E3, page S29')
#ggsave("figs/LocalMoransMBE3_2.pdf", plot=last_plot(), width=5.67, height=3.3)
## 2.11 Global Moran's I statistics, Member bills (Table E2) ----
table.format.mb(Gmorans.mb, 'g')
## 2.12 Spatial lag results, Member bills (Table E3) ----
table.format.mb(spatlag.mb, 's')
## 2.13 Distribution of nonunanimous votes (Figure E4) ----
merged %>%
  filter(!is.na(pctyea), pctyea<1) %>% 
  ggplot(.) +
  geom_histogram(aes(pctyea), color='black',fill='white', binwidth = .1, boundary = 0) +
  scale_y_continuous(limits=c(0,1000)) +
  labs(title = 'DISTRIBUTION OF NONUNANIMOUS VOTES',
       subtitle = 'ref: Appendix Figure E4',
       x = 'Percent Yea', y = 'frequency')
#ggsave("figs/NonUnVotesE4.pdf", plot=last_plot(), width=5.5, height=4, units = 'in')
## 2.14 Lopsided voting on all measures (Table E4) ----
merged %>%
  mutate(yea9 = case_when(pctyea>=0.9 & pctyea<=1 ~ 1),
         yea8 = case_when(pctyea>=0.8 & pctyea<.9 ~ 1),
         yea7 = case_when(pctyea>=0.7 & pctyea<.8 ~ 1),
         yea6 = case_when(pctyea>=0.6 & pctyea<.7 ~ 1),
         yea5 = case_when(pctyea>=0.5 & pctyea<.6 ~ 1),
         yea4 = case_when(pctyea>0.4 & pctyea<=.5 ~ 1),
         yea3 = case_when(pctyea>0.3 & pctyea<=.4 ~ 1),
         yea2 = case_when(pctyea>0.2 & pctyea<=.3 ~ 1),
         yea1 = case_when(pctyea>0.1 & pctyea<=.2 ~ 1), 
         yea0 = case_when(pctyea>=0 & pctyea<=.1 ~ 1), 
         indic = case_when(yea4==1 | yea5==1 ~ 'count', 
                           is.na(yea4==1) | is.na(yea5==1) ~ 'other')) %>%
  drop_na(pctyea) %>% 
  count(indic,session) %>% 
  pivot_wider(names_from = c(indic), values_from = n) %>% 
  mutate_if(is.numeric, ~replace_na(., 0)) %>%
  mutate('Percentage\n40-60%' = format(round(100*(count/(count+other)),2),nsmall=2), 
         'All Votes' = count+other) %>% 
  select(., -c(other, count)) %>% arrange(session) %>% gt() %>% 
  tab_header(title = "LOPSIDED VOTING ON ALL MEASURES", subtitle = 'rep: Appendix table E4') %>%
  tab_style(locations = cells_title(groups = c('title', 'subtitle')),
            style = list(cell_text( size='x-large', align = 'left')))